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ABSTRACT 

We present basic properties of protostellar disks in the embedded phase of star formation (EPSF), 
which is difficult to probe observationahy using available observational facilities. We use numerical 
hydrodynamics simulations of cloud core collapse and focus on disks formed around stars in the 
0.03 — 1.0 Mq mass range. Our obtained disk masses scale near-linearly with the stellar mass. The 
mean and median disk masses in the Class and I phases (M^g™ = 0.12 Mq, M^^g = 0.09 Mq and 
M™Qj" = 0.18 Mq, M™qj = 0.15 Mq, respectively) are greater than those inferred from observations 
by (at least) a factor of 2-3. We demonstrate that this disagreement may (in part) be caused by the 
optically thick inner regions of protostellar disks, which do not contribute to millimeter dust flux. We 
find that disk masses and surface densities start to systematically exceed that of the minimum mass 
solar nebular for objects with stellar mass as low as = 0.05 — 0.1 Mq. Concurrently, disk radii 
start to grow beyond 100 AU, making gravitational fragmentation in the disk outer regions possible. 
Large disk masses, surface densities, and sizes suggest that giant planets may start forming as early as 
in the EPSF, either by means of core accretion (inner disk regions) or direct gravitational instability 
(outer disk regions), thus breaking a longstanding stereotype that the planet formation process begins 
in the Class II phase. 

Subject headings: circumstellar matter — planetary systems: protoplanetary disks — hydrodynamics 
— ISM: clouds — stars: formation 



1. INTRODUCTION 

Theoretical and numerical studies indicate that spe- 
cific physical conditions are needed for planets to form 
in circumstellar disks and understanding disk properties 
may therefore help to choose between competing planet 
formation scenarios. In the past decade, much effort has 
been made to survey circumstellar disks around low- 
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Furlan et al.l 120091 and many others). However, these 



studies have mostly focused on the late evolution stage 
when the disk, if favorably located and inclined, is 
accessible for observations. On the contrary, the early 
embedded stage, when the disk is shrouded by a natal 
cloud core, is difficult to probe with the modern obser- 
vational techniques and only a ha ndful of attempted 
studies have been d one so far (e.g. lETsner et all 120051 : 
IChiang et all 120081 : iJorgensen eFall 120091 ). In this 
context, a numerical investigation into the properties 
of protostellar disks in the embedded phase of star 
formation (hereafter, EPSF) is of considerable interest. 

While there are plenty of theoretical and nu- 
merical studies addressing stability, fragmentation, 
and mass transport in cir c umste llar dis ks (e.g. 
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as masses, sizes, typical temperatures and densities have 
received considerably less attention owing to the (tech- 
nical) difficulty with forming circumstellar disks self- 
consistently in numerical hydrodynamics simulations. 
Indeed, numerical modeling of isolated disks cannot pro- 
vide reliable information on, for example, disk masses 
and sizes because they depend on the initial masses and 
rotation rates of parent cloud cores. 

On the other hand, numerical simulations that do 
form disks self-consistently often focus on the radial gas 
surface density and temperature profiles (|Lin fc PringlS 
1990t iHueso fc GuillotI [2005t iRice ct al. 2010; Vorobvo3 



2010al. accretion properties (e.g.|Dull emond et al. 200€ 



Vorobvov fc Ba su 2009) , and time evolution o f disk and 



stellar masses ()Nakamoto fc Nakagawal |1994| ) . In ad- 
dition, they often adopt many simplifying assumptions 
such as the barotropic equation of state, disk axisym- 
metry, etc. Three-dimensional numerical hydrodynamics 
simulations often suffer from a small number of consid- 
ered models due to aii enormous computa tional load (e.g. 
iMachida et~al]|2010l: iKratter et al.ll2010D . 

In this paper, we perform a comprehensive numerical 
analysis of embedded disks, which, for the first time, in- 
cludes obtaining the typical disk masses, radii, temper- 
atures, and densities in the Class and I phases of star 
formation for a wide spectrum of initial cloud core masses 
and angular momenta. We utilize two-dimensional nu- 
merical hydrodynamics simulations in the thin-disk ap- 
proximation with an accurate treatment of disk thermo- 
dynamics. This allows us to model the formation and 
long-term evolution of protostellar disks and derive use- 
ful statistical relations between various disk properties 
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and stellar masses and also make suggestions about the 
giant planet formation perspective in the EPSF. This is 
the third paper in a series focusing on the properties of 
circumstellar disks around (sub-)solar mass stars in the 
embedded phase of star formation. Two previous papers 
were mainly dedicat ed to the stability, fragm entation, 
accretion properties (jVorobvov fc Basrill2010bD, and also 
to the radial structure of protostellar disks (jVorobvovl 
l2010bD . 

The paper is organized as follows. In Section [5] we 
provide a brief description of our numerical model. Sec- 
tion [3] summarizes the initial parameters of cloud cores. 
Section |4] provides details on the adopted classification 
scheme for young stellar objects and also on the pro- 
cedure for distinguishing between the disk and infalling 
envelope. The main results are presented in Section [5] 
and are summarized in Section [T] The model caveats are 
discussed in Section [6] 

2. MODEL DESCRIPTION 

The m ain concepts of ou r appro ach are explained in 
detail in IVorobvov k, Basul ()2010bD and are briefly re- 
viewed below. We make use of numerical hydrodynamics 
simulations in the thin-disk approximation to compute 
the gravitational collapse of rotating, gravitationally un- 
stable cloud cores. This approximation is an excellent 
means to calculate the evolution for many orbital peri- 
ods and many model parameters. We start our numerical 
integration in the pre-stellar phase, which is character- 
ized by a collapsing starless cloud core, continue into 
the embedded phase of star formation, during which a 
star, disk, and envelope are formed, and terminate our 
simulations in the T Tauri phase, when most of the enve- 
lope has accreted onto the forming star/disk system. In 
the EPSF, the disk occupies the innermost region of our 
numerical grid, while the larger outer part of the grid 
is taken up by the infalling envelope, the latter being 
the remnant of the parent cloud core. This ensures that 
the protostellar disk is not isolated in the EPSF but is 
exposed to intense mass loading from the envelope. In 
addition, the mass accretion rate onto the disk Menv is 
not a free parameter of the model but is self-consistently 
determined by the gas dynamics in the envelope. 

We introduce a "sink cell" at Tsc — 6 AU and impose 
a free outflow condition on both boundaries so that the 
matter is allowed to flow out of the computational do- 
main but is prevented from flowing in. We monitor the 
gas density in the sink cell and when its value exceeds a 
critical value for the transition from isothermal to adi- 
abatic evolution 10^^ cm~'^), we introduce a central 
point-mass star. In the subsequent evolution, 90% of the 
gas that crosses the inner boundary is assumed to land 
onto the central star plus the inner axisymmetric disk at 
r < 6 AU. The other 10% of the accreted gas is assumed 
to be carried away with protostellar jets. 

Main physical processes that are taken into account 
in our modeling include stellar irradiation, background 
irradiation with temperature Tbg = 10 K, viscous and 
shock heating, radiative cooling from the disk surface, 
and also disk self-gravity. The corresponding equations 
of mass, momentum, and energy transport are 

By 

— = -V,.(St;,), (1) 



— (E-yp) + [V • {^vp ® = - VpT' + S + (2) 

+ (v-n)p, 

^ -I- Vp • [ev^) = -T'lVp • t;p) - A + F + [Wv)^^, : Hpp, , 

(3) 

where subscripts p and p' refer to the planar components 
(r, 0) in polar coordinates, E is the mass surface density, 

e is the internal energy per surface area, V = Pdh 
is the vertically integrated form of the gas pressure P, 
h is the radially and azimuthally varying vertical scale 
height determined in each computational cell using an 
assumption of local hydrostatic equilibrium, Vp — Vrf" + 

v^4> is the velocity in the disk plane, g„ = grV + g4,4> 
is the gravitational acceleration in the disk plane, and 

Vp — rd/dr + 4)r~^d/d4> is the gradient along the planar 
coordinates of the disk. 

Viscosity enters the basic equations via the viscous 
stress tensor 11. We parameterize the magnitude of 
kinematic viscosity v using a modified form of the a- 
prescription 

= aCshTair), (4) 

where Cg — is the square of effective sound speed 

calculated at each time step from the model's known V 
and S. The function J"„(r) = 27r-nan-i [(rd/r)i°] is a 
modification to the usual a-prescription that guarantees 
that the turbulent viscosity operates only in the disk and 
quickly reduces to zero beyond the disk radius . In this 
paper, we use a spatially and temporally uniform a, with 
its value set to 5 x 10~^. 

Radiative cooling from the disk surface is determined 
using the diffusion approximation of the vertical radia- 
tion trans port in a one-zone model of the vertical disk 
structure pohnson fc Gammiell2003l ) 

where a is the Stefan-Boltzmann constant, T is the mid- 
plane temperature of gas, and J-c — 2-|-20tan~^(r)/(37r) 
is a function that secures a correct transition between 
the cooling function (from both surfaces of the disk) in 
the optically thick regime Athick = IGctT'^/St and the 
optically thin one Ath in = '^'^'^'^J- We use frequency- 
integrated opacities of iBell fc LinI (|1994[ ). 
The heating function is expressed as 

F ^ F^a ,^ 2 , (6) 

L -\- T 

where Tirr is the irradiation temperature at the disk sur- 
face determined by the stellar and background black- 
body irradiation as 

^ irr ^ bg ^ : I ' J 

where Tbg is the uniform background temperature (in our 
model set to the initial temperature of the natal cloud 
core) and F{„{r) is the radiation flux (energy per unit 
time per unit surface area) absorbed by the disk surface 
at radial distance r from the central star. The latter 
quantity is calculated as 

Fii.r(r) = Aii.,.-^COS7ii.r, (8) 
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where 7irr is the incidence angle of radiation arriving 
at the disk surface at radial distance r, Ai„ is a time- 
dependent factor that is introduced to account for the 
possible attenuation of stellar radiation by the envelope'^, 
and is the sum of the accretion luminosity L*,accr aris- 
ing from the gravitational energy of accreted gas and the 
photospheric luminosity L*,ph due to gravitational com- 
pression and deuterium burning in the star interior. 

Viscous heating operates in the disk interior and is 
calculated using the standard expression (Vv)pp' : Tlppi. 
Heating due to shock waves is taken into account via 
compressional heating V (Vp • Vp) and artificial viscosity. 
The vertically integrated gas pressure V and internal en- 
ergy per surface area e are related via the ideal equation 
of state V = (7 — 1) e, with the ratio of specific heats 
7 = 7/5. 

Equations (H])-® are solved using the method of finite 
differences with a time-explicit, operator-split solution 
procedure in polar coordinates (r, (f>) on a numerical grid 
with 512 X 512 grid zones. Advection is treated using the 
van Leer interpolation scheme. The update of the inter- 
nal energy per surface area e due to cooling A and heating 
r is done implicitly using the Newton-Raphson method 
of root finding, complemented by the bisection method 
where the Newton-Raphson iterations fail to converge. 
The radial points are logarithmically spaced. The in- 
nermost grid point is located at the position of the sink 
cell = 6 AU, and the size of the first adjacent cell 
varies in the 0.07-0.1 AU range depending on the cloud 
core size. This corresponds to the radial resolution of 
A r-=1. 1-1.6 AU at 100 AU More details can be found 
in lVorobvov fc Basul ()2010b[ ). 

3. INITIAL CONDITIONS 

Initially isothermal (Tjnit = Tbg — 10 K) cloud cores 
have surface densities S and angular velocities fl typical 
for a collapsing, a xisymmetric, magnetically supercritical 
core (IBasulll997D 




where VIq is the central angular velocity and tq is the ra- 
dius of central near-constant-density plateau defined as 
ro — VAcI / (ttCEo) . With this choice of ro, equation ([9]) 
at large radii r ^ tq leads to the gas volume density 
distribution p = Acg/(27rGr^), if the latter is integrated 
in the vertical direction assuming a local vertical hydro- 
static equilibrium, i.e., p = T,/(2h) and h = Cg/{TrGY}). 
This means that our initial gas surface density configu- 
ration can be considered to have a factor of A positive 
density enhancement compared to that of the singular 
isothermal sphere psis — Cs/(27rGr^). Throughout the 
paper, we use A = 1.2. 

We present results from 4 model sets, the parameters 
of which are summarized in Table [1] Every model set 

^ The effect of this attenuation factor on the disk evolution is 
found to be insignificant, since A^^^ quickly approaches unity with 
time. For a typical run, Ai^^ 0.75 at t = 0.1 Myr after the 
formation of the central star and A^^^ Ri 0.95 at t = 0.3 Myr. 



TABLE 1 
Model parameters 



Model set 


/3 




Uo 


ro 




N 


1 


2.7 X 10" 


-S 


0.5-2.0 


960-3940 


0.4-1.8 


6 


2 


5.6 X 10 


-2 


0.7-6.0 


445-3770 


0.2-1.7 


8 


3 


1.3 X 10- 


-2 


1.2-12 


340-3430 


0.15-1.5 


8 


4 


2.3 X 10- 


-2 


2.1-29 


190-2570 


0.085-1.2 


7 



Note. — All masses are in Mq, distances in AU, and angular 
velocities in km pc~^. 

has a distinct ratio /3 — -Erot / 1 -Bgrav | of the rotational to 
gravitational energy calculated as 

l^out Tout 

Ey-ot = 2Tr J rOcSrdr, iJgrav = — 27r J rgr'Srdr. 

Here, Oc = fl^r and g-r are the centrifugal and gravita- 
tional accelerations, respectively, and Tout is the core's 
outer radius. The adopted v a lues o f /3 lie within the lim- 
its inferred by iCaselli et al.l ()2002D for dense molecular 
cloud cores, 10~^-7 x 10~^. In addition, every model set 
is characterized by a distinct ratio rout/?'o = 6 in order 
to generate gravitationally unstable truncated cores of 
similar form. As a result, individual cloud cores within 
each set of models have equal /3 and rout / ro but distinct 
masses Md, outer radii Tout, and central angular veloci- 
ties Qq. 

The actual procedure for generating a specific cloud 
core with a given value of /3 is as follows. First, we 
choose the outer cloud core radius Tout aud find tq from 
the condition rout/'^o = 6. Then, we find the central 
surface density Sq from the relation tq = V^c'^/{ttGT,q) 
and determine the resulting cloud core mass Md from 
Equation ([9]). Finally, the central angular velocity flo is 
found from the condition /3 = 0.9iloro/Cs. In total, we 
have simulated numerically the time evolution of 29 cores 
spanning a range of initial masses between 0.085 Mq and 
1.8 Mq. The adopted initial core m ass functio n is sim - 
ilar to that presented in figure 1 of iVorobvovl (l2010al) . 
The effect of different initial E and f2 configurations are 
considered in Section [01 

4. CLASSIFICATION SCHEME AND DISK-TO-ENVELOPE 
TRANSITION BOUNDARY 

Modern classification schemes of young stellar objects 
(YSOs) are designed to distinguish between main physi- 
cal phases of the early stellar evolution and the spectral 
energy distribution is often used to relate a YSO to a 
particular class (see lEvans et al.l[2009L for a thorough re- 
view). In our case, it is more con yenient to use the c lassi- 
fication breakdown suggested bv I Andre et al.l ()1993() and 
based on the mass remaining in the envelope 

Class Monv > 0.5Mci, 

Class I O.lMci < Afenv < 0.5Mci, (12) 

Class II Monv < O.lMoi- 

According to this scheme (hereafter, AWTB scheme), 
transition between Class and Class I objects occurs 
when the envelope mass Monv decreases to half of the 
initial cloud core mass Md. The Class II phase ensues 
by the time when the infalling envelope nearly clears and 
its total mass drops below 10% of the initial cloud core 
mass. Of course, we acknowledge that these numbers 
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Fig. 1. — Radial profiles of the azimuthally averaged gas surface 
density (top panels), radial velocity (middle panels), and angular 
velocity (bottom panels) at t = 0.07 Myr (left column) and t = 
0.11 Myr (right column) after the formation of the central star. The 
vertical dashed lines mark the position of a tentative surface density 
threshold for disk to envelope transition S = 0.1 g cm~^ (top 
panels) and the location of the disk outer edge as signalized by the 
characteristic deviation of the radial velocity from that expected 
for a free fall motion (dotted lines in middle panels). The dotted 
line in the bottom panels illustrates a Keplerian rotation. 

are somewhat arbitrary and the use of other classifica- 
tion diagnostics may introduce a systematic bias in our 
results. 

In order to use the AWTB scheme, we need to know 
which part of our numerical grid is occupied by the disk 
and which part belongs to the infalling envelope at any 
time instance of the evolution (note that due to the use 
of the sink cell we know the stellar mass). This is not 
trivial and we describe the adopted method in detail be- 
low. In order to distinguish between the infalling enve- 
lope and burgeoning disk, we make use of the gas sur- 
face density distribution in the computational grid. We 
set a threshold density for transition between the disk 
and envelope at I]d2c = 0.1 g cm~^, motivated by the 
observational fact that protoplanetary disks often de- 
cline exponentially w ith radius at S < 0.1 g cm~^ (e.g. 
[Andrews et al. I [20091) . Numerical hydrodynamics simu- 
lations of disk structure and evolution seem to confirm 
this phenomenon (Vorobyov 2010b)- 

Of course, using only I]d2e is not sufficient and some 
physically motivated quantity is to be invoked in order 
to minimize possible uncertainties. Therefore, we also 
make use of the gas radial velocity v^. The envelope is 
freely falling onto the disk and we determine the radial 
position where this free-fall motion terminates near the 
disk outer edge. 

To illustrate our method, we consider a typical model 
core with Md = 1.23 A/©, Sq = 4.5 x 10"^ g cm'^, tq = 
2750 AU, Q.0 = 1.25 km s^^ pc^^, and rout = 0.08 pc. 
Figure [1] shows (from top to bottom) the azimuthally 
averaged radial profiles of E, v^, and Q. dX t — 0.07 Myr 
(left column) and t = 0.11 Myr (right column) after the 



formation of the central star. The vertical dashed lines 
indicate an adopted transitional gas surface density of 
Sd2e = 0.1 g cm~^ (top panels) and the radial position 
where the infall motion of gas terminates at the disk 
outer edge (middle panel). The latter effect is mani- 
fested by a sharp deviation of the radial velocity from 
the typical free-fall profile cx r^'^-^ plotted by the dot- 
ted line. It is seen that the radial locations marked by 
the vertical dashed lines differ insignificantly from each 
other. It is worth noting that using Keplerian rotation 
arguments may not be that fruitful as one might expect. 
Indeed, the bottom panel demonstrates that a transition 
from Keplerian rotation, shown by the dotted line and 
typical for the disk, to sub-Keplerian one, typical for the 
infalling envelope, proceeds smoothly with a fairly wide 
transitional region. The use of this criterion may intro- 
duce significant uncertainties in the definition of the disk 
radius. 

To summarize, our procedure consists of the following 
steps. We march from the outer computational bound- 
ary toward the inner one and identify the radial loca- 
tion at which the azimuthally averaged gas surface den- 
sity is greater than Ed2e- Simultaneously, we determine 
the location where the envelope material hits the disk 
outer edge. We take the minimum of these two values 
as the location of the disk outer edge. To the best of 
our knowledge, this method for discriminating between 
the disk and envelope is found to be most accurate for 
our numerical model (thin disk). One may argue that 
we need not to use Sd2c at all, given that the other cri- 
terion is physically motivated and sufficient. We, how- 
ever, honestly believe that the exponential decline of the 
gas surface densi ty at large radii, as observed in cir- 
cumstellar disks ([Andrews et al.ll2009[ ). is not acciden- 
tal but physically motivated by, e.g., photoevaporation 
(HoUcnbach at al. 2000) , tidal stripping of the disk outer 
regions due to cl ose encounters between the members of 
a stella r cluster ([Bate fc Bonnelj [20031) . or viscous dis- 
persal (jVorobvov fc Basu[[2009b[ ). It is therefore impor- 
tant to impose some cut-off value in E to take these ef- 
fects into account. We acknowledge that our choice of 
Sd2c = 0.1 g cm~^ may somewhat affect the resulting 
disk masses and radii and we consider the effect of vary- 
ing Ed2c in more detail in Section [6.21 

Finally, we note that some of our model cores with 
high enough mass and angular momentum form bi- 
nary/multiple systems via disk fragmentation. When- 
ever this happens, we have to determine if the companion 
is still embedded in the disk of the primary or it has de- 
tached from the natal disk and assembled a circumstellar 
disk of its own. We do this by scanning the azimuthally 
averaged gas surface density and velocity profiles in our 
computational domain. The companion's position is usu- 
ally marked by a sharp peak in E and local flow of gas 
pointing to the companion. We determine if these two 
conditions are satisfied. In addition, we postulate that 
the companion completely detaches from its natal disk if 
a gap develops between the natal disk and the compan- 
ion with azimuthally averaged E < 0.01 g cm~^. Accord- 
ing to our tests, this procedure produces reliable results. 
When the companion forms and detaches from the natal 
disk, the mass of the latter drops accordingly. This pro- 
cess can be seen in the upper-right panel of Figure |8] in 
Section [6?T] 
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Fig. 2. — Disk masses versus stellar masses for the Class objects 
(top) and Class I objects (bottom) in model set 1 (open circles), 
model set 2 (filled squares), model set 3 (open triangles up), and 
model set 4 (filled triangles down). The solid lines show the least- 
squares best fits to the model data. 

5. INTEGRATED DISK PROPERTIES 

5.1. Disk masses 

A systematic numerical study of d isk masses around 
stars of (sub-)solar mass was done bv IVorobvovl (|2009a[ ) 
in the context of polytropic disks. In this paper, we have 
implemented a more accurate prescription for the disk 
thermal physics and expanded our parameter space by in- 
cluding cloud cores with higher initial rotational energies. 
Another important improvement over our previous work 
is the realization that both disk and stellar masses in- 
crease significantly during the EPSF. Therefore, present- 
ing the Md-M, relation in terms of the time-averaged 
quantities, as was done in ,Vorobyovi ( :2009af ), may not be 
adequate. In this study, we have chosen six typical disk 
and stellar masses based on the ratio of the envelope to 
cloud core mass, Monv/A^ci- This ratio is useful because 
it gradually decreases with time and can be used as a 
tracer of the early evolution of a protostar. In particu- 
lar, disk and stellar masses at Mcnv/-^^ci=0.85, 0.7, and 
0.51 are meant to represent typical ones in the begin- 
ning, in the midway, and at the end of the Class phase, 
respectively. For the typical disk and stellar masses in 
the Class I phase, we have chosen three values based on 
-Wcnv/-^^ci=0.49, 0.3, and 0.1. Of course, these six disk 
and stellar masses cannot represent the whole possible 
spectrum of masses that can be observationally detected 
in the embedded phase of star formation. Nevertheless, 
they can help to reduce a possible bias toward a par- 
ticular evolution stage and can help to make our disk 
mass versus stellar mass relation observationally mean- 
ingful. Finally, we have taken into account the mass of 



gas contained in the sink cell. The inner inflow compu- 
tational boundary is located at rsc=6 AU. The gas that 
flows through this boundary is assumed to land onto the 
central star and inner circumstellar disk with radius Tsc ■ 
The mass partition between these two objects can affect 
the obtained stellar and disk masses and is estimated 
by extrapolating the azimuthally averaged gas surface 
density profile into the inner 6 AU. Unfortunately, this 
procedure was not originally incorporated into the code, 
but we perform postprocessing our data to calculate the 
mass contained in the inner 6 AU. On average, the es- 
timated mass partition between the star and the inner 
disk is 10:1, with the effect of lowering the stellar mass 
by about 10% and increasing the total disk mass by the 
corresponding amount. 

Figure [2] presents our model disk masses versus stel- 
lar masses in the Class phase (top) and Class I phase 
(bottom) for four model sets, the parameters of which 
are summarized in Tabic [T] In particular, model set 1 is 
plotted by open circles, model set 2 — by filled squares, 
model set 3 — by open triangles up, and model set 4 — 
by filled triangles down. Each symbol of same shape 
within a given set of models represents an individual ob- 
ject formed from a core of distinct mass, rotation rate, 
and outer radius. Our modeling covers a wide range of 
central object masses, staring from sub-stellar objects 
with ~ 0.03 Mq and ending with solar-type stars. 

Several interesting features can be seen in Figure [5J 
First, there is an obvious correlation between the disk 
and stellar masses in the EPSF. The least-squares best 
fits to the Classes and I data (solid lines) yield the 
following relations 

Md^co = (0.73l°:J^)M:;°fO-"^ (13) 

Md,CI=(0.65^ro^)Ml■o^o^ (14) 

where the subscripts CO and CI refer to the Class and 
Class I phases, respectively, and disk and stellar masses 
are in solar masses. It is evident that the correlations 
are slightly super-linear and do not depend significantly 
on the particular phase. This result is in some disso- 
nance with what was found ea rlier in the context of poly- 
tropic disks (| Vorobvovll2009al ) . where a steepening of the 
Md-M* correlation was reported along the Classes O-II 
evolutionary sequence. We believe that this steepening 
was partly caused by time averaging of the correspond- 
ing disk and stellar masses over the duration of each 
phase, which effectively reduced the range of disk and 
stellar masses available for constructing the best fits. It 
may also be partly caused by a smaller parameter space 
in our earlier study, which span a narrower range of 
/3 = (1.2-3.4) X 10-3. 

Observations of protostellar disks in the EPSF are ex- 
tremely difficult due to the obscuration of light by sur- 
rounding envelopes and also due to uncertainties in dust 
opacities, optical depths, and disk structure. Disk masses 
inferred using complicated disk plus envelope models 
that employ Monte-Carlo radiative transfer codes and 
fit various circumstellar dust distributions to the mea- 
sured spectral energy distribution and (sub-)millimeter 
continuu m emission exist only for a handful of ob- 
jects (e.g.lBrown et a l. 2000; Andrew s fc Williamsll20"05t 
lEisner et al.ll2005l: |j0rgensen et al..,2009[ ). On the other 
hand, the T Tauri phase lacks notable envelopes and disk 
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masses in this phase are inferred for a number of objects. 
Yet, there is conflicting evidence as to the relation be- 
tween dis k and stehar masse s in the T Tauri phase. For 
instance, iNatta et all (|200lD found a marginal correla- 
tion between and Af*, albeit with a substantial dis- 
pers ion. On the o t her ha nd,, Andrew s fc Williams (.2005,) 
and lEisner et al.l (|2008D claimed no correlation. If we 
assume by extrapolation the existence of a similar (to 
Equations ([T3l) and (fT4|)) correlation between Md and A/» 
in the T Tauri phase, then our theoretical predictions are 
more in line with what was found by Natta et al, (,2001j). 

The numerically obtained near-linear correlation be- 
tween disk and stellar masses can be broken if a large 
population of objects occupies both the upper left and 
lower right portions of the M^-M^, phase space in Fig- 
ure [2] While it is feasible that there exists a fraction of 
stars with disk-to-star mass ratios low enough to popu- 
late the lower-right corner of Figure [2] (but see discus- 
sion in Section 15. 2p . it is unlikely that an equal frac- 
tion of sub-stellar objects with disks considerably more 
massive than the star could populate the upper-left cor- 
ner of Figure [2] The exist ence of such systems finds 
little observational (see e.g. [Andrews &: Williams! l2005f ) 
and theoretical s upport (see e.g. Kratter et al.ll2010D . As 
demonstrated by IVorobvovl ()2010a|) . svstems with equal 
disk and stellar masses are short-lived and quickly evolve 
into binary /multiple systems with disk masses consider- 
ably smaller than those of the stars. We conclude the 
our M^-Mf, correlation may weaken somewhat if more 
model cores with low values of /3 are considered (such 
cores are expected to form low-mass disks due to small 
centrifugal radii), but it cannot vanish completely. 

In accordance with our earlier results, we see a wide 
scatter in the derived disk masses for stars of equal mass, 
in particular for those in the intermediate mass range 
0.05-0.5 Mq. This scatter is caused by the fact that 
cloud cores with higher rotation rates (as defined by j3 
in our models) tend to form more massive disks. This 
implies that the distribution in initial conditions of cloud 
cores results in a scatter of disk masses. 

Table [2] presents our model mean, median, and max- 
imum disk masses (M^'^'^", M^'^", and Af™'''', respec- 
tively) calculated from the data shown in Figure [2j Our 



derived values are M™^^ 



0.12 Mq and M, 



mdn 
d,CO 



0.09 Mp 



o 



in the Class phase and also A^^ci" 



0.18 Mq and M^^^ = 0.15 Mq in the Class I phase. 
The existing estimates seem to yield lower mean disk 
masses, al beit with a consid erable scatter. For instance, 
[Andrews fc Williams! (|2005D reported M'^^f = 0.03 Mq 
for Class I sourc es in the Taurus- A uriga star forma- 
tion region, while !Brown et all (j2000! ) found Af^^g^^" = 
0.01 Mq for Class objects in the Perseus and Ser- 
pens molecular clouds. A recent sub-millimeter survey 
of low-mass protostar s in the Class and I phases by 
iJorgensen et al.l (j2009f ) revealed a somewhat higher mean 
disk mass of 0.05 Mq, yet a factor of 2-3 lower than 
our mean values. The same tendency of observationally 
inferred disk masses being smaller than those obtained 
from numerical simulations seems to extend to the later 
phases of stellar evolution (jVorobvov!!2009al ) . 

Given the difficulties with the observational determi- 
nation of disk masses in the EPSF and a low number 
statistics of embedded sources, we believe that drawing 




0.2 0.4 

Time (Myr) 

Fig. 3. — Time evolution of disk masses in model 1 (top), model 2 
(middle), and model 3 (bottom). In particular, the solid lines 
represent the total disk mass M^, while the dashed lines show 
the mass M^(t < 1) of those disk regions that are characterized 
by the frequency-integrated optical depth to the midplane r < 1. 
The dash-dotted lines are the ratio Md(T < 1)/Md. 



any firm conclusions from the comparison of our numeri- 
cal results with the observations at this stage may be pre- 
mature until more powerful observational facilities (like 
ALMA) come into work and more observational data be- 
come available. Nevertheless, we would like point out 
one possibility that may explain why the observation- 
ally inferred disk masses seem to be systematically lower 
than our model estimates. We believe that this discrep- 
ancy is at least partly caused by high opacity of the inner 
disk regions, which may cont ain a sizeable fraction of the 
total mass budget ()Rice et a l. 2010; Vorobyov 2010a). 
This possibility has a lso been put forward recently by 
iGreaves fc Ricg (l2010l) . 

To estimate the magnitude of the mass deficit due to 
high disk opacity, we calculate the disk mass that is char- 
acterized by the optical depth to the midplane r < 1 and 
compare the resulting value with the total disk mass. We 
consider three typical cores with distinct initial masses: 
model 1 is characterized by Mc\ = 0.2 Mq, while mod- 
els 2 and 3 have M^i = 0.85 Mq and M^i = 1.7 Mq, re- 
spectively. The ratio /3 is identical for these models and 
is set to 5.6 x 10"'^. Figure [3] shows the total disk mass 
(Md, solid lines), the mass of disk regions with t < 1 
{Md{T < 1), dashed lines), and the ratio A/d(T < 1)/Afd 
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TABLE 2 
Averaged disk properties 



Phase 


/\ /mean 


1\ -fmdn 


]\ /max 


^mcan 


^mdn 


^max 




~,mdn 


Class 


0.12 


0.09 


0.45 


0.71 


0.67 


1.30 


230 


140 


Class I 


0.18 


0.15 


0.53 


0.66 


0.66 


1.18 


480 


290 



Note. — Disk masses are in Mq and radii arc in AU. 



(dash-dotted lines) as a function of time since the on- 
set of gravitational collapse. In particular, the top panel 
shows the results for model 1 (top), while the middle 
and bottom panel present the data for model 2 and 3, 
respectively. The vertical dotted lines mark the onset of 
the Class I phase (left) and Class II phase (right). The 
value of M^{t < 1) is used as a proxy to what can be ex- 
pected from the measurements of disk masses in (partly) 
optically thick protostellar disks. 

As anticipated, Md(r < 1) is always smaller than 
and the difference is particularly large in the embedded 
phase but tends to diminish in the Class II phase. The 
rate of convergence between Md{T < 1) and is how- 
ever different in models with distinct core masses. The 
ratio Mfi{T < 1)/Md approaches unity much slower in 
higher- Afci models, indicating that massive disks remain 
partly optically thick for a longer time. On a more qual- 
itative side, Md(T < 1) in the Class phase may be 
a factor of 5-10 smaller than the total disk mass. In 
the Class I phase, the difference is less impressive (usu- 
ally a factor of 2-3) but may become much stronger 
for low-mass disks. We note that we have excluded a 
possible contribution of the sink cell to the total disk 
mass budget due to a difficulty with calculating self- 
consistently the optical depth there. However, disk re- 
gions at r < = 6 AU are usually optically thick 
and this fact should further magnify the effect. We con- 
clude that, due to high opacity, the observationally in- 
ferred disk masses may be seriously underestimated in 
the Class phase and, to a lesser extent, in the Class I 
phase. 

Finally, we point out that both modern theories of gi- 
ant planet formation, those of core accretion and direct 
gravitational instability, require that the gas surface den- 
sity in circumstellar disks be several time s greater than 
that of the minimum mass solar n ebular (jPollack et alj 
[19961 IBossI [19981 llda fc Lin|[200l . Taking 350 AU as 
a characteristic disk radius (see Table ^ and adopting 
the M MSN gas surface density profile from lHavashi et ahl 
(|1985D . we obtain « 0.03 Mq for the mass of the MMSN. 
Figure [2] shows that most of our models have disks with 
masses in excess of 0.03 Mq (horizontal dotted line) in 
both the Class and Class I phases of star formation. 
The least-squares fits indicate that disk masses start to 
systematically exceed that of the MMSN for objects with 
stellar mass as low as A/* w 0.05 Mq. 

At a first glance, the abundance of objects with disk 
mass greater than that of the MMSN may suggest that 
giant planets may start forming as early as in the Class 
and I pha ses of stellar evolution, a supposition also ad- 
vocated by iGreaves fc Ricd ()2010[ ) based on the statis- 
tics of observationally inferred disk masses in different 
stages of star formation. This may certainly be so for 
the core accretion mechanism. However, the direct grav- 
itational instability scenario encounters serious difficul- 



ties in the EPSF. Although it is true that protostellar 
disks in the EPSF are more prone to gravitational in- 
stability and fragmentation than in any other stage of 
their evolution, it is also true that the conditions for 
gia nt planet survival are least fav orable in the EPSF. 
As IVorobvov fc Basul ()2006l . l2010bD have demonstrated, 
gravitational interaction of protoplanetary embryos with 
natal spiral arms results in rapid inward migration of the 
embryos to the inner few AU. The future prospects for 
the embryos are unclear due to the specifics of modeling 
(sink cell at a few AU), but judging from fast migra- 
tion timescales (a few orbital periods) we believe that 
most of them will be tidally destroyed and absorbed 
by the central star, leading to a luminosity outburst 
similar to that of the FU-Orionis-type or EX-Lupi-type 
stars. However, if dust sedimentation and H2 dissocia- 
tion timescales are comparable to that of the migration 
one (for instance, in massive enough embryos at wide 
enough orbits), then some of the embryos may survive 
this migration, loose their upper atmospheres and form 
either metal-reach giants or even earth-type planets on 
orbit s of order a few AU (iBolev et al.l 120101 : iNavakshinI 
120101 ICha fc NavakshinI |2010D. Alternatively, embryos 
that form in the very late EPSF are less exposed to rapid 
radial migration; they may open a gap and settle on wide 
orbits of order 100 AU (jVorobvov fc Basull2010al ). How- 
ever, the expected frequency of such systems is 10% at 
best. 

5.2. Disk-to-star mass ratios 

The disk-to-star mass ratio ^ = M^/M^ is another im- 
portant diagnostic that can give us some useful infor- 
mation about stability, fragmentation, and mass trans- 
port properties in circumstellar disks. Observationally 
inferred ^ vary in wide limits, with most of th e val- 
ues l y ing in the 0.01-0.05 range (e.g. [B cckwit h et al.l 
flggO; iMannings fc SargentI [2OOOI: [Andrews fc Williamsl 
2005). As to a possible dependence of disk-to-star 
mass ratios on stellar masses, this issue is unsettled 
and some authors claim no correlation (e.g. iNatta et al.l 
|2001|) . while others find a negative correlatioii. i.e. , 
highe r £ for low er M, (e.g. IMannings fc Sargentj I2OOOI : 
Andr ews fc Wil liams 20Q1). 

Figure 0] presents our model disk-to-star mass ratios 
derived from the data of Figure [5] for the Class (top) 
and Class I (bottom) objects. The meaning of the sym- 
bols is the same as in Figure [2J The mean disk-to-star 
mass ratios in the Class and I phases are ^co'*" = 0.71 
and fc^" ~ 0.66, respectively. These values are con- 
siderably greater than those inferred from observations 
owing to higher disk masses derived in our modeling. In 
addition, measurements of disk masses are mostly done 
for the Class II objects, for which we could expect lower 
values of ^ due to accretion of matter from the disk onto 
the star, disk viscous dispersal, and photoevaporation. 

Another interesting feature seen in Figure |4] is the ap- 
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parent lack of objects with ^ < 0.2. The matter is 
that our earhest measurements of disk masses are at 
Mcnv/Mc\ = 0.85. By this time instance, disks have al- 
ready gained some non-negligible mass owing to ineffi- 
cient inward mass transport in the very early phases of 
disk formation (when gravitational instability is underde- 
veloped and viscous torques alone cannot cope with mass 
loading from the infalling envelope). As a result, the disk 
quickly grows in mass until gravitational instability and 
fragmen tation set in and decelerate the process of mass 
growth ()Vorobvovll2009al l2010biV The net effect is that 
systems with ^ < 0.2 are expected to be statistically rare. 
In addition, observational estimates of ^ at the very early 
phases of star formation when Menv/Afci > 0.85 are very 
challenging. 

The least-squares best fits shown in Figure 2] by the 
solid lines yield the following relations between disk-to- 
star mass ratios in the Class and Class I phases (^co and 
^ci, respectively) and the corresponding stellar masses 
(in Mq) 



eco = (0.72 ± 0.1) Af°;°f"-°^ 



^ci = (0.64 ±0.04) AC^i 



0.01±0.04 



(15) 
(16) 



It is evident that ^ demonstrates little correlation with 
M^, in the EPSF. The maximum disk to star mass ra- 
tios found in our modeling in the Classes and I phases 
are ^co = 1-3 and fci — 1-18, res pectively. Th i s is in 
agreement with a recent study by iKratter et al.l ()2010l ) 
an d a factor of 3 gr eater than was previously reported 
bv IVorobvovl ()2009a|) in the context of polytropic disks. 

It is interesting that the maximum ^ is attained by 
objects with intermediate stellar masses M* ss 0.2 — 
0.3 rather than by the upper-mass stars in our 
sample (M* ss 1.0 Af,) as might have been intuitively 
expected. The matter is that the mass transport in 
disks around low- mass stars < 0.2 Mq is domi- 
nated by viscous torques rathe r than by gravitational 
ones (jVorobvov fc Basul l2009b[ ). As the stellar mass 
grows, so does the mass of the parent core and the 
viscous torques fail to transport a growing amount of 
infalling core material through the disk and onto the 
star. This causes ^ to increase in the 0.03-0.2 Mq mass 
regime. For stars with Af* > 0.2 — 0.3 Af*, the situation 
changes qualitatively — mass transport in their disks be- 
comes dominated by gravitational torques and, in even 
to a greater exten t, my quick inw ard migration of form- 
ing fragments t Vorobvov fc Basu 2010b). Both mech- 
anisms appear to be considerably more efficient mass 
transport agents than viscous torques and the corre- 
sponding disk-to-star mass ratios start to decline. An- 
other factor contributing to the decline in ^ is the forma- 
tion of binary /mu ltiple systems in massive enough disks 
(lVorobvovll2010bh . 

5.3. Disk radii 

The size of a circumstellar disk is an important charac- 
teristics containing indirect information about the disk 
stability properties. Numerous numerical and theoreti- 
cal studies indicate that circumstellar disks of larger size 
are more prone to the development of gravitational in- 
stability and fragmentation, in part due to an increased 
disk mass and in part due to faster radiative cooling and 
lesser stellar irradiation/ viscous heating at large radii. 
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Fig. 4. — Disk-to-star mass ratios versus stellar masses for the 
Class objects (top) and Class I objects (bottom) in model set 1 
(open circles), model set 2 (filled squares), model set 3 (open tri- 
angles up), and model set 4 (filled triangles down). The solid lines 
show the least-squares best fits to the model data. 

Indeed, the stellar irradiation fiux declines with radius 
as h/r^ (X r"^ ''^ for a flared disk with the typical ratio 
of the scale height to radius h/r (x r°'^^, which makes 
extended disks colder (for a fixed stellar luminosity) and 
hence more prone to gravitational instability. At the 
same time, the ratio tcooi/idyn of the characteristic cool- 
ing time icooi = e/A (where e is the internal energy per 
surface area and A is the cooling rate from the disk sur- 
face) to the dynamical time idyn = 27rfi~^ declines with 
radius as in the optically thick regime (for E cx r~^/^, 
T (X r~^/^, and spatially independent dust opacity) and 
becomes r-independent in the optically thin case, again 
suggesting that the outer disk regions are more suscepti- 
ble to gravitational instability. 

Figure [5] presents the time-averaged disk outer radii 
(rd) as a function of the time-averaged stellar masses 
(Af*) in the Class phase (top) and Class I phase (bot- 
tom). The meaning of the symbols is the same as in 
Figures [5] and [H We perform time averaging over the 
duration of the corresponding evolutionary phase in or- 
der to smooth out large rad ial pulsations see n in most 
disks during the EPSF (see IVorobvovl l2010bl for more 
detail). In order to inform the reader on the magnitude 
of these pulsations, we plot vertical bars representing the 
minimum and maximum disk radii found in each model. 
To avoid some ambiguity in the definition of the mini- 
mum disk radius, we calculate this quantity at the time 
when the envelope mass drops to 0.85% that of the initial 
cloud core mass, i.e., when A^cnv/Afci = 0.85. 

It is seen that, for a given stellar mass, disk radii 
may vary by as much as two orders of magnitude. The 
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mean and median disk radii in the Class phase are 
^d!co" = 230 AU and rJJ^^'J = 140 AU, respectively. The 
corresponding values in the Class I phase are r™^^" = 
480 AU and r^'^j = 290 AU. The characteristic radii in 
the Class I phase increase by a factor of 2 as compared 
to those in the Class phase, reflecting the ongoing disk 
growth and expansion due to infall of the envelope ma- 
terial during the EPSF. 

Figure [5] reveals a mild trend of protostellar disks to 
grow in size with increasing stellar mass (and disk mass 
due to the near- linear Afd~M, relation). The least- 
squares best fit to the model data yields the following 
relations between the time-averaged disk radii (in AU) 
and time-averaged stellar masses (in M©) 

(rd,co) = (450;?^0) ^M^^^^)0.,,±o.2,^ (^7) 

{r,,ci) = {850tlfo) {M*.cir''^'-'°, (18) 

where indices CO and CI refer to the Class and Class 
I phases, respectively. These sub-linear correlations can 
be understood if we assume that protostellar disks at- 
tain a self-regulatory state in which the disk surface 
density does not increase notably with increasing stel- 
lar mass. In this case, the mass of such disks will grow 
due to an increase in the disk size. This in turn implies 
rd oc M°-5 oc (for ^ = Md/M* oc M*, see FigureS]). 

A somewhat steeper dependence found in our numerical 
simulations is explained by the fact that in reality the 
disk mass grows due to an increase in both the disk size 
and gas surface density, though the former process domi- 
nates the latter in sufficiently massive disks forrned from 
cloud cores with mass Mci > 0.9 Mq (lVorobvovll20"l0bl ). 

Observational estimates of disk sizes in the EPSF en- 
counter the same difficulties as in the case of disk masses 
and only a few object s have been studied so far. For 
instance, lEnoch et al.l (|2009[ ) reported a massive disk 
around Serpens FIRS 1, a well-kn own Class s ource , 
with radius of order 300-500 AU. iHogerheiidd (|2001l ) 
found an even larger disk around L1489 IRS, presum- 
ably an object in transition between the embedded and 
T Tauri phases, with radius of order 2000 AU. More ob- 
servations of embedded sources are needed to confirm 
that large disks are typical for the embedded phase of 
star formation. 

The statistics is wider in the case of Class II disks. Ob- 
servational estimates of disk sizes for 1.0-Myr-old proto- 
planetary disks in the Orion star-forming region illumi- 
nated by the UV radiat ion of massive stars were done by 
IVicente fc AlvesI ()2005[ ). They found typical disk radii 
around a large sample of late-G to late-M stars in the 
Trapezium cluster to be in the 50-200 AU range, with 
a median value of 70 AU for protoplanetary disks and 
130 AU for a subset of silhouette disks. These estimates 
are conside rably smaller than our m odel values. On the 
other hand, iMann fc Willianisl (|2009l ) reported disk radii 
for two massive protoplanetary disks situated beyond the 
Trapezium cluster in Orion to be of order 300 AU, sug- 
gesting that strong UV radiation from OB stars may 
have evaporated part of the nearby disks and truncated 
their radii. This view is supported by the estimates of 
disk radii in the Taurus and Ophiuchus star forming re- 
gions, less extreme than tha t of the Orion, performed by 
[Andrews fc Williamsl ()2007D . They found disk radu to 
lie in a wide range between 50 AU and 1000 AU, with a 
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Fig. 5. — Time-averaged disk radii versus time-averaged stellar 
masses for the Class objects (top) and Class I objects (bottom) 
in model set 1 (open circles), model set 2 (filled squares), model 
set 3 (open triangles up), and model set 4 (filled triangles down). 
Vertical bars represent the minimum and maximum disk radii in 
each model. The solid lines show the least-squares best fits to the 
model data, while the dashed lines mark a fiducial critical radius 
for disk fragmentation via gravitational instability. 

median value of 200 AU. This value is in better agree- 
ment with our median values, which are r™QQ=140 AU 

for the Class disks and r™^'f=290 AU for the Class I 
disks. We stress, however, that the observationally in- 
ferred disk radii really correspond to a sensitivity limit 
to the surface brightness profile of the continuum emis- 
sion and real disks may be somewhat larger. 

Many theoretical and numerical studies indicate that 
disk fragmentation is unlikely in the inner 50-100 AU 
due to insufficient cooling and elevated viscous and 
stellar irradiation heatin g ()Rafikovl [20091: iClarkd [200l 
IVorobvov fc Basull2010b[ ). In the context of planet for- 
mation, this means that disks must be greater than 
100 AU in radius in order to form giant planets via 
disk fragmentation. Figure [5] reveals a substantial frac- 
tion of objects with time-averaged disk radii exceeding 
100 AU (as marked by the horizontal dashed lines), in 
the Class I phase in particular. These objects are formed 
from prestellar cores of sufficiently high initial mass and 
angular momentum. The least-squares fits indicate that 
disk radii start to exceed 100 AU for objects with stellar 
mass as low as 0.05-0.07 Mq. 

5.4. Characteristic gas surface densities and midplane 
temperatures 

In this section, we provide typical gas surface densi- 
ties and temperatures of protostellar disks in the em- 
bedded phase of their evolution. Due to the use of the 
sink cell, we can calculate typical values only at rela- 



10 



Vorobyov 



stellar mass (Mr.-)) 



Stellar mass (M,. 



Stellar mass (M,;;,) 



Stellar mass (M,. ) 



0.02 0.05 0.1 0.2 0.5 




Stellar mass {M,.^) 



Stellar mass (M,. 



Fig. 6. — Typical disk surface densities at 40 AU (left column) 
and 100 AU (right column) in the Class phase (top row) and 
Class I phase (bottom row). The solid lines provide the least- 
squares fits to the model data, while horizontal dotted lines mark 
the gas surface densities in the MMSN at 40 AU and 100 AU. 

tively large radii. This information is of particular inter- 
est for the giant planet formation mechanism via direct 
gravitational instability, which has been predicted and 
dem onstrated to operate in the disk outer regions only 
("e.g. iDodson-Robinson et al . 2009; Rafikov 2009; Clarke! 
l2009HVorobvov fc Basull2oTOa : Bolev et al.. ,2010). 

Figure El presents gas surface densities at 40 AU (S40, 
left column) and 100 AU (Eiooj right column) from the 
star as a function of stellar mass for the same four 
sets of model cores as in Figure [2] In particular, the top 
and bottom panels correspond to the data in the Class 
and Class I phases, respectively. As in Figure [5J the 
gas surface densities and stellar masses are calculated at 
the beginning, in the midway, and at the end of each 
evolutionary phase. Because we are interested in disk 
surface densities, we have excluded model cores that fail 
to form large enough disks at the evolutionary times of 
interest (about 10% of the total number of models). 

The least squares best fits (solid lines) yield the follow- 
ing relations 

(19) 
(20) 
(21) 
(22) 



^40, CO 



80+1° 



0.75±0.1 
CO ' 



^40,01 

^100, CO ~ 55- 
Eioo,ci = 10±2MO;^5±o-\ 



-3.5 ^"*,c/ 

-12 ,,r0.6±0.1 

-10^"*, CO ' 



where gas surface densities and stellar masses are in 
g cm~^ and solar masses. It is seen that the correlations 
between S and are sub- linear, more massive disks are 
generally denser than their low-mass counterparts, and 
E declines with radius as found in many theoretical and 
numerical studies. 

As discussed in Section 15.11 the gas surface density 
should be several times greater than that of the minimum 
mass solar nebular (MMSN) for the gas giants to form 
either through core accretion or gravitational instability. 
The horizontal dashed lines in Figure[6]mark the gas sur- 
face densities of the MMSN at 40 AU and 100 AU as ca l- 
culated from the radial profile of iHavashi et al.l ()1985[ ): 
SiviMSNfe cm-2] ^ 1700(r/l AU)'^-^. The gas surface 
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Fig. 7. — Typical gas midplane temperatures at 40 AU (left 
column) and 100 AU (right column) in the Class phase (top row) 
and Class I phase (bottom row). The solid lines provide the least- 
squares fits to the model data. 

density at 100 AU is most relevant to giant planet for- 
mation via direct gravitational instability, which is sup- 
posed to operate at distances of order 100 AU and be- 
yond. The least-squares fits indicate that gas surface 
densities at 100 AU start to systematically exceed that 
of the MMSN for objects with stellar mass as low as 
0.07-0.09 Mq. 

On the other hand, the core accretion scenario is be- 
lieved to operate in the inner 10 AU. We cannot provide 
reliable data for these inner regions due to the specifics 
of our modeling. However, we note that as one goes from 
100 AU to 40 AU the corresponding gas surface densi- 
ties start to exceed that of the MMSN for objects of 
even smaller stellar mass M*=0.03-0.05 Mq (in contrast 
to 0.07-0.09 Mq for Sioo)- If we assume by extrapola- 
tion that this trend continues to even smaller radii, than 
planet formation via core accretion can also occur in ob- 
jects with stellar mass as low as 0.03-0.05 Mq. 

Finally, in Figure [7] we provide typical gas midplane 
temperatures at 40 AU {T^q^ left column) and 100 AU 
(Tioo, right column) for the same models as in Figure El 
The top/bottom rows show data for the Class 0/Class I 
phases, respectively. The least squares best fits (solid 
lines) yield the following relations 



T4o,co = 90±8M°■™^ (23) 
^4o,CI = 50±5M°■5±o■"^ (24) 

^loo,co = 55±4Mr™^ (25) 

rioo,ci = 30±5M°-™5. (26) 

It is seen that massive disks (i.e., those around stars of 
greater mass) are generally hotter than their low-mass 
counterparts. On the other hand, the correlations be- 
tween T and Af* are notably weaker than those between 
S and M*, indicating that efficient cooling operates at 
radii of order 50-100 AU. There seem to be a shght de- 
crease in the exponents along the Class 0-Class I evolu- 
tionary sequence. One interesting feature of Figure [7] is 
a notable jump in the gas midplane temperature seen at 
around A/* — 0.05 — 0.1 A/q, which represents a transi- 
tion from optically thin to optically thick disks. 
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Fig. 8. — Top. Disk masses (dashed lines), stellar masses (thin 
solid lines), envelope plus companions masses (dash-dotted lines), 
and disk-to-star mass ratios (thick solid lines) versus time passed 
since the formation of the central star in model 4 (left) and model 5 
(right). Bottom. Disk radii versus time in model 4 (left) and 
model 5 (right). The vertical dotted lines mark the onset of the 
Class I phase (left) and Class II phase (right). 

6. MODEL CAVEATS 

6.1. Initial conditions 

Our model cloud cores have distinct masses (correlated 
with variations in Sg and external pressure) and rota- 
tional energies (correlated with fto). At the same time, 
the shape of the radial gas surface density and angular 
velocity profiles remains fixed (Equations © and PU)) ). 
Although the E oc profile is appropriate for super- 
critical cores formed via ambipolar difi^usion and, per- 
haps, for a wider fa mily of cores formed via slow gravi- 
tational contraction (jPapp fc Basu 2009), it is not guar- 
anteed, however, that this profile is universal. In par- 
ticular, the current model assumes each star is roughly 
the Jeans mass in its natal environment and, as a con- 
sequence, more massive stars come from a lower pres- 
sure environment. However, the assumed independence 
of the shape of prestellar cores from the external pres- 
sure may break, for massive star formation in particu- 
lar^. The same arguments apply to the adopted angu- 
lar velocity pr ofile. As an a l ternat ive, we take the ap- 
proach of Bos s fc HartmannI (|2001| ) and explore an ini- 
tially isothermal, self-gravitating, sheetlike cloud with 
volume density depending only on distance from the mid- 
plane p{z) = p{0)sech^ {z / h) . Such cylinder- like con- 
figurations have constant surface densities and are also 
solid-body rotators. In the following, we consider two 
specific model prestellar cores that are characterized by 
equal masses Md = 1.23 Mq and ratios /? = 9 x 10~^ 
but distinct initial shapes. In particular, model 4 has 
nonuniform S and f2 described by Equations ([9]) and 
(fTO]) . with Eo = 4.5 X 10~2 g cm"2^ m = 2750 AU, and 
f^o = 1.25 km s^^ pc^^, while model 5 has spatially uni- 
form E = 1.3 X 10^-^ g cm~^ and 17 = 1.25 km s~^ pc~^. 
In both models, rout — 0.08 pc. These two models 
are meant to represent two limiting cases for the ini- 
tial shapes of prestellar cores in low-mass star forming 
regions. 

The top panels in Figure |8] present stellar masses 

We note, however, that our stellar mass spectrum is limited 
from above by stars with mass roughly equal to that of the Sun. 
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TABLE 3 

TiME-AVERAGED MODEL CHARACTERISTICS IN CLASS PHASE 

Model (gco) (Md^co) {M,,co) (A/cnv,co) (''d.co) 

4 0.71 0.13 0.18 0.9 325 

5 0.87 0.17 0.19 084 453 

Note. — All masses are in Mq and disk radii — in AU. 

TABLE 4 

TiME-AVERAGED MODEL CHARACTERISTICS IN CLASS I PHASE 

Model (gci) {Mi,ci) (Ah^ci) (A/cnv,Cl) {rj'^ 

4 0.63 0.34 0.55 0.27 1065 

5 0.64 0.33 0.52 0.29 1060 

Note. — All masses are in Mq and disk radii — in AU. 

(thin solid line), disk masses Md (dashed lines), envelope 
plus companion masses M^nv (dash-dotted lines), and 
disk-to-star mass ratios ^ (thick solid lines in model 4 
(left) and model 5 (right) as a function of time passed 
since the formation of the central star. The bottom pan- 
els show the disk radius in model 4 (left) and model 5 
(right). The vertical dotted lines mark the onset of 
Class I and Class II phases of star formation. The two 
models show distinct time evolution, as expected for 
gravitationally unstable protostellar disks formed from 
non-identical prestellar cores, but when time averaged 
values are considered it turns out that notable differ- 
ences take place only in the Class phase. Tables [3] and 
[4] present main model characteristics time-averaged over 
the duration of Class and Class I phases, respectively. 
From Table [3] it is evident that the nonuniform model 5 
has greater ^po, ^'^d,C0j Af*,co, and rd,co but smaller 
-^^onv.co than the uniform model 4, indicating that pro- 
tostellar cores with spatially uniform profiles of E and 
n sustain higher infall rates onto the disk and deplete 
their material faster than cores with spatially declining 
profiles. At the same time. Table |4] reveals little dif- 
ference between the two models in the Class I phase. 
We conclude that, if the whole spectrum of initial condi- 
tions were taken into account, we might expect to obtain 
somewhat higher disk masses, stellar masses, disk-to-star 
mass ratios, and disk radii in the Class phase but the 
corresponding values in the Class I phase would remain 
largely unaffected. 

6.2. Discriminating between the disk and injalling 
envelope 

Throughout the paper, we have adopted a character- 
istic transitional density between the disk and envelope 
of Ed2e = 0.1 g cm^^. Here, we illustrate the effect 
of a factor of two smaller Ed2e = 0.05 g cm~^. Top 
panel in Figure [9] shows the disk mass versus time for 
Ed2c = 0.1 g cm~^ (sohd line) and Ed2c = 0.05 g cm~^ 
(dashed line), while the bottom panel presents disk radii 
for the corresponding transitional densities. It is evident 
that the disk mass is little affected by a factor of 2 drop 
in Ed2c, owing to an exponentially declining gas surface 
density in the disk outer regions. The disk radii are some- 
what more sensitive to Ed2c, yet showing only a factor 
of 1.2 increase on average. Figure [S] convincingly demon- 
strates that our model results are not critically sensitive 
to the choice of Ed2e- 
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Fig. 9. — Top. Disk mass versus time for Sd2c = 0.1 g cm~^ 
(solid line) and 2^20 = 0.05 g cm~^ (dashed line). Bottom. Disk 
radii for the corresponding vales of Sd2G- 

6.3. Jet efficiency and limitations of the thin-disk 
approximation 

In this paper, we have taken the simplest approach and 
assumed that 10% of the mass passing through the sink 
cell is subsequently ejected with protostellar jets. This 
number has some motivation since observations of the 
variation of jet diameters with distance from their driving 
sources have been consistent with m odels giving the jet 
efficiency of > 3% (|Rav et al.ll2"007l) . At the same time, 
theoretical and numerical studies suggest a higher jet ef- 
ficiency, of up to 70% depending on the degree of cor e 
flattening (e.g. lShu et al.lll994l : lMatzner fc McKeil2000D . 
due to a three-dimensional nature of the jet phenomenon. 
In this case, not only the wind material is lost but also 
the envelope material is entrained and evacuated, driving 
massive molecular outflows. It is however very difficult 
to self-consistently take into account the removal of the 
envelope material by jets in our two-dimensional simu- 
lations. We postpone the solution of this problem to a 
later time, when an accurate vertical disk structure is 
implemented in the code, and note that a higher jet effi- 
ciency (than currently adopted) may result in lower final 
stellar and disk masses. 

Another issue associated with the thin-disk model is 
that the envelope material is forced to land all at once 
onto the disk outer edge, while in three dimensions it 
may cover a larger area above and below the disk sur- 
face. However, numerical simulations of gas trajectories 
in flaring disks (our model disks have the aspect ratio 
h/r cx r°'^^) indicate that the bulk of the envelope ma- 
terial (up to 80%) still lands onto the disk outer edge 
(jVisser et al.l 120091 ). The remaining fraction falls onto 
the disk surface, triggering the streaming instability due 
to the vertical sheer in the gas velocity. This instabil- 



ity acts to further destabilize the disk and increase the 
mass and angular rnomen tum transport in protostellar 
disks ()Harsono et al.|[201ClD . an effect missing in the cur- 
rent two-dimensional disk model. The net result is some 
decrease in the disk mass and an equivalent increase in 
the stellar mass. 

7. CONCLUSIONS 

Using numerical hydrodynamics simulations in the 
thin-disk approximation, we model the formation and 
evolution of protostellar disks in the embedded phase of 
star formation (EPSF) when the disks are exposed to in- 
tense mass loading from natal cloud cores. A wide spec- 
trum of initial core masses and ratios /3 of the rotational 
to gravitational energy is considered, which allows us to 
construct statistical relations between disk masses Md, 
disk outer radii r^, and stellar masses M*. Typical gas 
surface densities and midplane temperatures at 40 AU 
and 100 AU are derived and analyzed as a function of 
stellar mass. Our findings can be summarized as follows. 

• Our numerical modeling reveals that most of the 
Class and I objects with stellar masses > 
0.05 — 01 Mq are characterized by disk masses 
and gas surface densities greater than that of 
the minimum m ass solar n ebula r, the latter de- 
fined as in Hava shi et al.l ()1985f ). At the same 
time, disk radii for these objects start to grow 
beyond 100 AU, making gravitational fragmenta- 
tion possible. These results suggest that giant 
planet formation may commence as early as in the 
Cl ass phase, subs t anti ating an earlier conclusion 
of Grea ves fc Ricd (|2010 ) based on the statistics 
of observationally inferred disk masses in different 
stages of star formation. 

• Disk masses in the EPSF show a near-linear cor- 
relation with stellar masses. The correlation may 
weaken somewhat if more models with low (3 are 
considered but is unlikely to break completely. 

• A wide scatter of disk masses for a given stellar 
mass is likely caused by the distribution in the ini- 
tial conditions of cloud cores and, in particular, by 
a wide spectrum of initial rotation rates. 

• Our mean and median disk masses in the Class 
and I phases (M|;g^" = 0.12 Mq, M^^gg = 0.09 Mq 

and M™gf = 0.18 Mq, M^gj" = 0.15 Mq, respec- 
tively) are greater than those inferred from obser- 
vations by (at least) a factor of 2-3. We show that 
this difference may be caused by high opacity in 
the inner disk regions. 

• Disk-to-star mass ratios are found to lie mostly in 
the 0.2 < £, < 1.0 range, with the mean and median 
values ^ w 0.65 — 0.7 in the Class and I phases. 
Objects with ^ < 0.2 are rare due to fast accumu- 
lation of disk mass in the early Class phase, while 
objects with ^ > 1.0 quickly evolve into binary or 
multiple systems, thus effectively reducing the disk 
mass. The least-squares fitting reveals little corre- 
lation between ^ and M* in the EPSF. 
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• Our modeling predicts a sub-linear (close to a 
square root) scaling between time-averaged disk 
outer radii and stellar masses. Disks grow in size 
in the EPSF due to mass loading from collapsing 
cores. Our mean and median values in the Class 



and I phases (rg^g^o" 
and rJJ^g'j" 



230 AU, rj}^^^ = 140 AU 



480 AU, r'l'^l 



290 AU, respectively) 
suggest that protostellar disks are large enough for 
giant planets to form via disk fragmentation. 

• We find sub-linear correlations between gas surface 
densities at 40 AU and 100 AU and stellar masses 
and also a square-root scaling between disk mid- 
plane temperatures and stellar masses, implying 
that massive disks are denser and hotter than their 
low-mass counterparts. 

Although disk conditions seem to be favorable for gi- 
ant planets to start forming as early as in the EPSF, 



we want to stress that the direct gravitational insta- 
bility scenario may encounter serious difficulties caused 
by rapid inward migrations of the fo rming protoplane- 
tary e mbryos jV orobyov & Basu 2006:. [20"l0bl : iNavakshinI 
120101 : ICha fc Na vakshin .2010) . Accor ding to recent nu- 
merical hydrodynamics simulations bv lVorobyov fc Basul 
(|2010al) . the fraction of systems with survived embryos 
amounts to 10% at best. 
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